#!/bin/bash
set -e

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------

###--------------- argument may be changed ---------------###



# genome_name=b37.fa
echo Using $genome_name as the reference genome 2>>$log 1>&2
# genome_assembly=b37
echo genome assembly is: $genome_assembly 2>>$log 1>&2
# dbsnp_version=138
echo dbsnp version is: $dbsnp_version 2>>$log 1>&2

# data_thread_num=10
# cpu_thread_num=4

# java_memory=-Xmx16g
# java_memory=16g
echo java memory: $java_memory 2>>$log 1>&2

# snpeff_db_version=GRCh37.75 
echo $snpeff_db_version 2>>$log 1>&2

# header=\
# @HD\tVN:1.4\tGO:none\tSO:coordinate

# read_group=\
# @RG\\tID:${sample_name}\\tPL:ILLUMINA\\tSM:$sample_name
# For gatk, Each read group must contain the platform (PL) and sample (SM) tags.
# For the platform value, we currently support 454, LS454, Illumina, Solid, ABI_Solid, and CG (all case-insensitive).
# Each read in the file must be associated with exactly one read group.

vcf_path=${data_path}/vcf/gatk/

# there is no space between "=" and "\"
# ref_genome=\
# ${data_path}/ref/b37/$genome_name

echo ref genome is: $ref_genome



###--------------- end argument may be changed ---------------###





# java -jar snpEff.jar download -v GRCh37.71


# 经 gatk 和 snpeff 处理后，anno FILTER 报错，所以先运行
# ${script_path}anno_fp_0424.sh 1.snp.vcf 1.snp.fp.vcf 2>>$log
# ${script_path}anno_fp_0424.sh 1.indel.vcf 1.indel.fp.vcf 2>>$log


echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $gatk VariantFiltration 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.snp.filtered.vcf \
	--variant 1.snp.vcf \
	--filterExpression "FS > 60.0" \
	--filterName "FS_more_than_60" \
	--filterExpression "MQ < 40.0" \
	--filterName "MQ_less_than_40" \
	--filterExpression "HaplotypeScore > 200.0" \
	--filterName "HaplotypeScore_more_than_200" \
	2>>$log
	
echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $gatk VariantFiltration 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.indel.filtered.vcf \
	--variant 1.indel.vcf \
	--filterExpression "FS > 200.0" \
	--filterName "FS_more_than_200" \
	2>>$log



echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $gatk CombineVariants 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T CombineVariants \
	--variant 1.snp.filtered.vcf \
	--variant 1.indel.filtered.vcf \
	-o 1.snp.indel.filtered.vcf \
	-genotypeMergeOptions UNSORTED \
	2>>$log
	
# echo
# echo
# echo VEP annotation
# variant_effect_predictor.pl \
# --offline \
# --force_overwrite \
# --dir $data_path/ensembl/ \
# --fasta $ref_genome \
# --refseq \
# --hgvs \
# --pick \
# --format vcf \
# --vcf \
# --input_file 1.snp.indel.filtered.vcf \
# --output_file 1.filtered.vep.vcf



# Annotate using ENCODE's data:
# java -Xmx4g -jar snpEff.jar -v -interval db/encode/wgEncodeDukeDnase8988T.fdr01peaks.hg19.bb GRCh37.75 examples/test.1KG.vcf > test.1KG.ann_encode.vcf



echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift annotate id 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
annotate \
-id ${vcf_path}dbsnp_${dbsnp_version}.${genome_assembly}.vcf \
1.snp.indel.filtered.vcf \
> 1.snpsift.dbsnp.vcf \
2>>$log
# -tabix               : VCF database is tabix-indexed. Default: false



echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift annotate id 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
annotate \
-id ${data_path}/ncbi/clinvar/clinvar_20150106.vcf \
1.snpsift.dbsnp.vcf \
> 1.dbsnp.clinvar.vcf \
2>>$log




echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift annotate cosmic 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
annotate \
${data_path}/cosmic/CosmicCodingMuts.vcf \
1.dbsnp.clinvar.vcf \
> 1.id.cosmic.vcf \
2>>$log

echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift annotate varType 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
varType \
1.id.cosmic.vcf \
> 1.cosmic.var_type.vcf \
2>>$log


echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift annotate gwasCat 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
gwasCat \
-db ${data_path}/snpeff/gwascatalog.txt \
1.cosmic.var_type.vcf \
> 1.var_type.gwas.vcf \
2>>$log

:<<!

echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift dbnsfp 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
dbnsfp \
-db ${data_path}/dbnsfp/dbNSFP2.9.txt.gz \
-f Ensembl_transcriptid,Uniprot_acc,Interpro_domain,SIFT_score,Polyphen2_HVAR_score \
1.var_type.gwas.vcf \
> 1.gwas.dbnsfp.vcf \
2>>$log

!

# echo 2>>$log 1>&2
# echo 2>>$log 1>&2
# echo $snpeff ann snpeff2gatk 2>>$log 1>&2
# time \
# java -Xmx$java_memory -jar $snpeff \
	# -c ${snpeff_path}/snpEff.config \
	# -v \
	# -lof \
	# -noStats \
	# -o gatk \
	# $snpeff_db_version \
	# 1.gwas.dbnsfp.vcf \
	# > 1.dbnsfp.snpeff2gatk.vcf \
	# 2>>$log
	
	# -t 8 \ ??
	
	
	
echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpeff ann snpeff2gatk 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpeff \
	-c ${snpeff_path}/snpEff.config \
	-v \
	-lof \
	-noStats \
	$snpeff_db_version \
	1.var_type.gwas.vcf \
	> 1.dbnsfp.snpeff2gatk.vcf \
	2>>$log
	
echo
echo
echo VEP annotation
variant_effect_predictor.pl \
--offline \
--force_overwrite \
--dir $data_path/ensembl/ \
--fasta $ref_genome \
--refseq \
--hgvs \
--pick \
--format vcf \
--vcf \
--input_file 1.dbnsfp.snpeff2gatk.vcf \
--output_file 1.vep.vcf


:<<!
echo
echo
echo $gatk VariantFiltration
java -Xmx$java_memory -jar $gatk \
	-R $ref_genome \
	-T VariantFiltration \
	-o 1.gatk.filter.vcf \
	--variant 1.dbnsfp.snpeff2gatk.vcf \
	--filterExpression "QD < 2.0" \
	--filterName "gatk_QD_less_than_2" \
	--filterExpression "MQRankSum < -12.5" \
	--filterName "gatk_MQRankSum_less_than_-12.5" \
	--filterExpression "ReadPosRankSum < -8.0" \
	--filterName "gatk_ReadPosRankSum_less_than_-8" \
	--filterExpression "ReadPosRankSum < -20.0" \
	--filterName "gatk_indel_ReadPosRankSum_less_than_-20" \
	--filterExpression "InbreedingCoeff < -0.8" \
	--filterName "gatk_indel_InbreedingCoeff_less_than_-0.8"
	
!





# cat 1.dbnsfp.snpeff2gatk.vcf \
# | awk -F', ' '/##INFO=<ID=\w+, /{print $1","$2","$3","$4}'

cat 1.vep.vcf \
| awk -F', ' '{if (/##INFO=<ID=\w+, /) { print $1","$2","$3","$4}else{print $0}}' \
> 1.header_modified.vcf



echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo $snpsift filter 2>>$log 1>&2
time \
java -Xmx$java_memory -jar $snpsift \
filter \
"(FILTER = 'PASS')" \
1.header_modified.vcf \
> 1.high.vcf \
2>>$log


# echo 2>>$log 1>&2
# echo 2>>$log 1>&2
# echo $snpsift filter 2>>$log 1>&2
# time \
# java -Xmx$java_memory -jar $snpsift \
# filter \
# "(FILTER != 'PASS')" \
# 1.dbnsfp.snpeff2gatk.vcf \
# > 1.low.vcf \
# 2>>$log



$bgzip -c 1.header_modified.vcf > 1.header_modified.vcf.gz
$tabix 1.header_modified.vcf.gz

echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo bcftools view for filter FILTER 2>>$log 1>&2
time \
$bcftools view \
-f "FP_VARIANTCALLER,FP_PRIMERDIMER" \
1.header_modified.vcf.gz \
| bgzip -c \
> 1.fp.vcf.gz \
2>>$log

gunzip -c 1.fp.vcf.gz > 1.fp.vcf

$tabix 1.fp.vcf.gz

$bgzip -c 1.high.vcf > 1.high.vcf.gz
$tabix 1.high.vcf.gz
$bcftools concat -aD 1.high.vcf.gz 1.fp.vcf.gz |bgzip -c > 1.high_fp.vcf.gz
$tabix 1.high_fp.vcf.gz

if [ ! -d isec ]; then
	mkdir isec
fi
$bcftools isec -p ./isec 1.header_modified.vcf.gz 1.high_fp.vcf.gz

cp ./isec/0000.vcf 1.low.vcf



# echo 2>>$log 1>&2
# echo 2>>$log 1>&2
# echo bcftools view for filter FILTER 2>>$log 1>&2
# time \
# $bcftools view \
# -f "FS_more_than_60,MQ_less_than_40,HaplotypeScore_more_than_200,FS_more_than_200" \
# 1.dbnsfp.snpeff2gatk.vcf \
# > 1.low.vcf \
# 2>>$log



# java -jar SnpSift.jar rmInfo test.snpeff.vcf EFF
